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Abstract 

The voltage trace of neuronal activities can follow nnultiple timescale dynamics that arise from correlated membrane 
conductances. Such processes can result in power-law behavior in which the membrane voltage cannot be characterized 
with a single time constant. The emergent effect of these membrane correlations is a non-Markovian process that can be 
modeled with a fractional derivative. A fractional derivative is a non-local process in which the value of the variable is 
determined by integrating a temporal weighted voltage trace, also called the memory trace. Here we developed and 
analyzed a fractional leaky integrate-and-fire model in which the exponent of the fractional derivative can vary from 0 to 1, 
with 1 representing the normal derivative. As the exponent of the fractional derivative decreases, the weights of the voltage 
trace increase. Thus, the value of the voltage is increasingly correlated with the trajectory of the voltage in the past. By 
varying only the fractional exponent, our model can reproduce upward and downward spike adaptations found 
experimentally in neocortical pyramidal cells and tectal neurons in vitro. The model also produces spikes with longer first- 
spike latency and high inter-spike variability with power-law distribution. We further analyze spike adaptation and the 
responses to noisy and oscillatory input. The fractional model generates reliable spike patterns in response to noisy input. 
Overall, the spiking activity of the fractional leaky integrate-and-fire model deviates from the spiking activity of the 
Markovian model and reflects the temporal accumulated intrinsic membrane dynamics that affect the response of the 
neuron to external stimulation. 
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Introduction 

The leaky integrator properties of a neuron are determined by 
the membrane resistance and capacitance which define a single 
time constant for the membrane voltage dynamics [1-3]. 
However, the voltage trace of real neurons can follow multiple 
timescale dynamics [4—6] that arise from the interaction of 
multiple active membrane conductances [7-11]. Such processes 
can result in power-law behavior in which the membrane voltage 
cannot be characterized with a single time constant [5,12-15]. 
Since power-law dynamics can span all the scales of interest of 
neuronal behavior [16-20], it is necessary to develop a framework 
to study such processes and their effect on the electrical and 
computational capacities of neurons. 

In the classical leaky integrate-and-fire model the temporal 
evolution of the voltage is local [21,22]. The value of the voltage at 
a given time depends only on the value of the voltage in the 
immediate previous time step. Such a process is called Markovian. 
However, coupling of active conductances does not allow the value 
of the voltage to be memoryless [1 1,17,23-26]. Instead, long time 
correlations affect the membrane voltage for hundreds of 
milliseconds. The emergent effect of these membrane correlations 
is a non-Markovian process that can be modeled with a fractional 
derivative [27-31]. A fractional derivative represents a non-local 
process [32-34] in which the value of the variable is determined by 



integrating a temporal weighted voltage trace, also called the 
memory trace. Although fractional derivatives and integrals are 
almost as old as traditional calculus [32,35,36], they have not been 
widely used due to limited computer power. In the fractional 
integrate-and-fire model the exponent of the fractional derivative 
goes from 0 to 1, with 1 representing the normal derivative. As the 
exponent of the fractional derivative decreases, the weights of the 
voltage trace increase. Thus, the value of the voltage is increasingly 
correlated with the trajectory of the voltage in the past. 

We developed and analyzed a fractional leaky integrate-and-fire 
model. The only parameters of the model are the conductance, 
capacitance, and the fractional exponent. By varying the fractional 
exponent our model reproduces the upward and downward spike- 
frequency adaptations found experimentally in pyramidal neurons 
[37,38], tectal neurons [39] and fast-spiking cells of layer IV in the 
rat barrel cortex [40]. Furthermore, the model replicates not only 
the adapting firing rate but also the long first-spike latency seen in 
pyramidal neurons in layer V [38]. The model also produces 
spikes with longer first-spike latency and high inter-spike 
variability with power-law distribution, which cannot be repro- 
duced by the classical integrate-and-fire model. We further 
analyze spike adaptation and the responses to noisy and oscillatory 
input. Overall, the spiking activity of the fractional integrate-and- 
fire model deviates from the spiking activity of the Markovian 
model and reflects the temporal accumulated intrinsic membrane 
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Author Summary 

Spike adaptation is a property of most neurons. When 
spike time adaptation occurs over multiple time scales, the 
dynamics can be described by a power-law. We study the 
computational properties of a leaky integrate-and-fire 
model with power-law adaptation. Instead of explicitly 
modeling the adaptation process by the contribution of 
slowly changing conductances, we use a fractional 
temporal derivative framework. The exponent of the 
fractional derivative represents the degree of adaptation 
of the membrane voltage, where 1 is the normal leaky 
integrator while values less than 1 produce increasing 
correlations in the voltage trace. The temporal correlation 
is interpreted as a memory trace that depends on the 
value of the fractional derivative. We identify the memory 
trace in the fractional model as the sum of the instanta- 
neous differentiation weighted by a function that depends 
on the fractional exponent, and it provides non-local 
information to the incoming stimulus. The spiking 
dynamics of the fractional leaky integrate-and-fire model 
show memory dependence that can result in downward or 
upward spike adaptation. Our model provides a framework 
for understanding how long-range membrane voltage 
correlations affect spiking dynamics and information 
integration in neurons. 

dynamics that afiect the response of the neuron to external 
stimulation. 

Results 

The objective of this project was to develop a fractional leaky 
integrate-and-fire model of neuronal activity to study spiking 
adaptation. For this purpose we developed a fractional differential 
model of the leaky integrator combined with a regular spiking 
generation mechanism. Complex multiple timescale neuronal 
systems can be studied using fractional or power-law dynamics; 
examples range from ion channel gating properties, to diffusion of 
intracellular signals in Purkinje and pyramidal cells, synaptic 
strength and firing rate adaptation [14,19,27,37,41-44]. 

The fractional leaky integrate-and-fire model 

We define the fractional leaky integrate-and-fire model as 



-gL(V-VL) + h, 



along with the fire-and-reset condition 



if V{t) = F,h, then spike at time t and V^^ 



(2) 



where V is the membrane potential, and a is the order (exponent) 
of the fractional derivative, with 0 < a < 1 . In the case of a = 1 , the 
fractional model is the same as the classical leaky integrate-and-fire 
model. When the membrane potential reaches a threshold (Kth), a 
spike is generated and V is reset to Kreset for a refractory period 
Tref. The passive membrane time constant is 
Parameter values are given in Table 1 (see Methods). For 

0 < a < 1 , the fractional derivative of the voltage (-— — in Eq . 1 ) 
can be defined with the Caputo [45] fractional derivative 



d^V _ 1 
liF ~ r(l-0£). 



0 



dr. 



(3) 



By numerically integrating the above fractional derivative (Eq. 
3) using the LI scheme [46], we approximate the fractional 
derivative of order a, where 0<o;< 1, 



d''V (dty 



dF r(2-a) 

'N-l 

Y,[^('k+i)-ntkmN-kf-^^-{N-\-kf-''^] 

k = 0 



(4) 



where F(/o)=Fo, and tk is the value of time such that 
tk = kdt. For all simulations, we use the time step dt = Q.\ ms. By 
combining the right sides of Eqs. 1 and 4, and solving for V at time 
In [V{tfj)) that depends on all past values of V (from K(?o) to 
V{tfi-\)), we obtain 



V{ti^)K(dtfY(2-(y) 



-gL(VitN-i)-VL) + h.] 

Cm 



+ n'iv-i) 



Y^[V{tk+i)-V{tkmN-kf-'''^-{N-\-kf-''\ 



(5) 



where we define the Markov term weighted by the gamma 
function as 



{dtfT{2-a) 



-gL{V{tM-l)-VL) + hn 



Cf}] 



+ V(t„_iX (6) 



and the voltage-memory trace as 

J2[^('''+o-nh)m-kf-''-^-(N-i-kf-''^]. (7) 



The voltage-memory trace (Eq. 7) can be further divided into 
the differentiation of past voltage {AV{tic)) (Eq. 8) weighted by a 
function Wi^(a,k) that depends on o; (Eq. 9): 

AV(tk)= V(tk+i)-V(tk), for /t = 0,1, 2, ■■■,N-2, and (8) 



WN{oi,k) = {N-kf-''^-(N-l-kf-''\ 
for A: = 0,l,2, •••,Af-2, 



(9) 



where k is the time counter for past events and the positive integer 
corresponds to the value of time tf/ (/at = Ndt) at which the 
voltage V(tf^) is integrated (Eq. 5). 

The voltage-memory trace contains information of aU the 
previous voltage activity of the neuron. Clearly, this is a 
computationally intensive problem due to the expanding matrix 
over time. We integrate this equation using our recently developed 
Fractional Integration Toolbox [47]. Since the fractional integra- 
tion in Eq. 5 needs at least two inputs, V(ti) is first integrated 
using the classical leaky integrate-and-fire model. 
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Table 1. Parameter values for the fractional leaky integrate-and-fire model. 





Parameter 


Value 


Description 


Cm 


0.5 nP 


membrane capacitance 




-70 mV 


leak reversal membrane potential 


v,^ 


-70 mV 


resting membrane potential 


Preset 


-70 mV 


reset membrane potential 




-70 mV 


initial membrane potential 




-50 mV 


threshold membrane potential 




25 nS 


leak conductance 


km 


3 nA 


injected (applied) current 




5 ms 


refractory period 



gL = 16.7 nS that results in t,„ = 30 ms. 
doi:1 0.1 371 /journal.pcbi.l 003526.t001 



i and 



The fractional leaky integrate-and-fire model shows spike 
adaptation 

In this section we show tlie spiking properties of the fractional 
leaky integrate-and-fire model and compare our results to 
experimental data, mainly from cortical pyramidal neurons. 

Spike adaptation depends on the value of the fractional 
order. We start by comparing the classical and fractional leaky 
integrate-and-fire models. The classical leaky integrator model is 
fijlly described by the cell membrane resistance and capacitance 
(Fig. lA left). As mentioned before we propose that the emergent 
contribution of slowly varying conductances results in voltage 
dynamics that can be modeled with a fractional leaky integrator 
(Fig. lA right). We first compared the voltage response of the 
model without the spiking mechanism, thus studying only the sub- 
threshold dynamics. In such conditions the voltage response of the 
fractional model to a step current using a = 1 is identical to the 
response of the classical model (Fig. IB left). However, as a 
decreases the voltage response shows adaptation that is not 
characterized by T„, (Fig. IB right). Instead, the voltage response 
increasingly shows power-law adaptation as the value of 0£ 
decreases (Fig. IC). Thus, the fractional derivative transforms 
the sub-threshold voltage response to a constant stimulus from 
exponential to a power-law. 

We characterized the spike adaptation of the model as a 
function of a by stimulating the model with a step current and 
measuring the response to the first spike and the properties of the 
inter-spike intervals (ISIs). Our results show that the sub-threshold 
voltage dynamic of the fractional model is reflected in the spiking 
activity of the neuron, with no adaptation and identical spiking 
activity as the classical model when a = 1 and increasing spike 
adaptation as a decreases (Fig. ID). Our model shows an inverted 
dependency of the latency of the first-spike to the value of oc 
(Fig. IE). The first-spike latency is common in the activity of 
neurons and has been suggested as a source of information for 
decision accuracy [48-50]. The ISI of the model in response to 
constant stimulation shows increasing adaptation as the value of 0( 
decreases. An analysis of the ISI responses shows that for values 
a<0.2 the slope of the ISI follows a power-law 
distribution(Fig. IF). Power law distribution of ISIs may maximize 
frring rate entropy [51]. It is important to emphasize that this 
variability arises from the intrinsic properties of the model since 
the input is constant. As a consequence our fractional order model 
produces spike trains with a wide range of inter-spike intervals, 
with power-law dynamics for smaller values of a (Fig. IG). 



The formulation of the fractional leaky integrate-and-fire model 
and the results in Fig. 1 suggest that the activity of the neuron 
depends on the integration of previous activity beyond the value of 
T„, . If the neuronal dynamics depend on the memory trace then 
the amplitude and time of a previous input would aflFect the spiking 
activity of the fractional model. To determine this point we 
injected two different amounts of a hyper-polarizing current for a 
fixed period followed by an identical depolarizing current (Fig. 2A). 
As is well known, the firing rate response of classical leaky 
integrate-and-fire model shows basically identical responses except 
for the time to first spike after the onset of the depolarization. 
However, the fractional model suggests a more complex response 
in which not only the time of the first spike but the instantaneous 
firing rate of the ensuing spiking response depends on the value of 
the hyper-polarizing current (Fig. 2B). Similar results are obtained 
when there is a fixed hyper-polarizing current applied for different 
amounts of time. In this case, the longer the fractional model 
integrates the hyper-depolarizing stimulus, the more it affects the 
delay to first spike and spike adaptation (Fig. 2C). 

Since the response of the fractional model is history dependent, 
the firing rate curve versus injected current might change 
depending on previous input. To test this possibility we injected 
either a hyper-polarizing (pre-stimulus low. Fig. 2D top, dashed 
lines) or depolarizing (pre-stimulus high. Fig. 2D top, solid lines) 
current to the model before injecting current for 6 s to calculate 
the firing rate. We calculated the mean firing rate of the spike train 
produced during the last 3 s of the stimulation. This quantification 
of the firing rate shows that as a decreases the threshold to 
generate spiking activity increases. However, the firing rate curves 
resulted in the same value independent of the previous stimulation 
(Fig. 2E, the dashed hues for the pre-stimulus low and solid lines 
for the pre-stimulus high overlap), thus showing that after 3 s the 
fractional model reaches steady state while showing a typical Type 
I dependency [52,53]. Although the steady state response of the 
system is independent of past activity, the initial instantaneous 
firing rate of the first inter-spike interval shows dependence on 
previous activity [54]. For o;< 1.0, the initial firing rate calculated 
from (pre-stimulus high) (Fig. 2F, solid) is always higher than the 
initial firing rate calculated from (pre-stimulus low) (Fig. 2F, 
dashed). Thus, this analysis shows that although the model is 
integrating all the voltage-memory trace the influence of this 
component decreases over time. 

Adaptation to periodic inputs. We decided to systemati- 
cally study the response properties of the fractional model to 



PLOS Computational Biology | www.ploscompbiol.org 



3 



March 2014 | Volume 10 | Issue 3 | el 003526 



Fractional Integrate-and-Fire Model 



Leaky inteqrate- 
and-fire 
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integrate-and-fire 
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Figure 1. Comparison between tlie classical and fractional leaky integrate-and-fire models. (A) Schematic circuit diagrams for the 
classical (left) and fractional order (right) leaky integrate-and-fire models. (B) Sub-threshold response in the classical (left) and fractional models (right). 
Both stimulated with /i„j = 0.3 nA. (C) The sub-threshold voltage response converges to a power-law function when a decreases. (D) While the 
classical model (left) generates regular spiking to a constant input, the fractional model (right) shows first spike latency and spike adaptation. Both 
models stimulated with /i,,, = 3 nA. (E) The first-spike latency produced by the fractional model becomes longer when a is smaller. (F) The inter-spike 
interval histogram as a function of a. The histogram shows power-law distribution as o!->0.2. (G) The inter-spike intervals decrease over time as a 
function of a. The color key in C applies to F and G. 
doi:10.1371/journal.pcbi.1003526.g001 



oscillatory input from the sub-threshold to the spiking regime. A 
common quantification method of the sub-threshold properties of 
a neuron is to use a frequency varying sinusoidal current, also 
known as a ZAP current (Fig. 3A top). As is well known [55,56], 
the amplitude of the sub-threshold voltage oscillations of the 
classical model decreases as the input current frequency increases. 
This is a consequence of the low-pass properties of the classical 
leaky integrator. Our model replicates this behavior for a = 1 . 
However, the voltage response decreases in amplitude as a 
decreases (Fig. 3A). We quantified the input-output relationship of 
the model by calculating the impedance magnitude in the 
frequency domain (|Z|) [57,58] (see Methods). As expected for 
the passive membrane model, the peak impedance is at the lowest 
frequency and there is no resonance for 0( = 1 . This property of the 
model is robust and not affected even for smaller values of oc. 
However, for smaller values of a. the impedance becomes almost 
constant across aU tested frequencies (Fig. 3B left). We also 
analyzed the phase shifts between the input current and the output 
voltage. The phases are always negative indicating lags in which 
the voltage oscillation follows the input current. However, the 
phase angle becomes less negative when a decreases indicating 
smaller lags which are practically the same across all frequencies 
(Fig. 3B right). Overall, our results show that lower values of a 
result in a low but homogeneous fdter of oscillatory inputs across 
all frequencies. 

In order to initially characterize the spiking response to an 
oscillatory input we used a sinusoidal current with constant 
frequency and amplitude (Fig. 3C). We chose parameters that 
generated supra-threshold spiking at the peak of the oscillation and 
no firing at the trough for a = 1 . As expected, as the value of oc 
decreases the response of the neuron is delayed. Interestingly, this 
results in a spiking response to the oscillatory input in which the 
number of spikes per cycle increases slowly. This behavior is due to 
the accumulation of the input in the memory trace of the fractional 



model (Fig. 3D). We then characterized the response of the neuron 

to supra-threshold oscillatory input (firing at any value of the input 

current) while varying the period of the stimulation. For this 

purpose we calculated the gain of the spiking activity of a neuron, 

which is defined as the ratio of the amplitude of the instantaneous 

firing rate to the amplitude of a sinusoidal input current with a 

fixed period length. Experimental results have shown that the 

firing rate of L2/3 pyramidal neurons behaves as a fractional 

differentiator [37] in which the firing rate follows the fractional 

/27i\" . f2nt o(7i\ 
derivative of the input. In this case, it is 1^1 ^''^ 1 ~2~ ) 

when the input is sin^^^-^ with period T. We simulated such an 

experiment in our model, with o£ = 0.15 (Fig. 3E). We then fitted 
the resulting instantaneous firing rate with a sine wave in which 
the free parameters were the amplitude and the phase. Our results 
show, as demonstrated experimentally, that the gain decreases as a 
function of the period following a power-law [37]. A least-square 
fit of these data points shows 

^gain — 0.15, which matches the order 
of the fractional leaky integrate-and-fire model used in the 
simulation (Fig. 3F). Although it has been reported that the phase 
lead is frequency independent [37], our results show that the phase 
lead decreases as the period of the sine wave increases (Fig. 3G), 
which is consistent with more recent reports [59]. The gain power- 
law behavior and phase lead response are preserved even in the 
presence of noise added to the sine wave input (not shown), again, 
replicating experimental results. This analysis shows that the 
fractional sub-threshold voltage behavior of the model is reflected 
in the spiking activity of the neuron. 

There is an increasing body of evidence showing that when a 
cortical neuron is stimulated with alternating current steps the 
spiking activity of the neuron shows adaptation dependent on the 
period of the stimulus [37-40]. Our model reproduces this 
dynamic, particularly when a<0.2 (Fig. 3). As done in the 
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Time after depolarization (s) I 

inj 



Figure 2. The mean and instantaneous firing rate responses of tKie fractional model to constant input. (A) Two levels of hyper-polarizing 
current followed by the same depolarizing current result in different spiking patterns. (B) The instantaneous firing rate against time for different 
conditions described in A. (C) The instantaneous firing rates against time for identical hyper-polarizing current (-3 nA) with different durations. (D-F) 
Comparison of mean and instantaneous firing rates. (D) Top: Applying a hyper-polarizing (pre-stimulus low, dashed lines) current before application 
of current steps to calculate firing rate responses. Bottom: as in the Top but applying a depolarizing (pre-stimulus high, solid lines) current. (E) The 
mean firing rates (mean FR) vs injected current show Type I response for both (pre-stimulus low) and (pre-stimulus high) current input paradigms 
described in D. The dashed and solid lines corresponding to pre-stimulus low and high, respectively, overlap. (F) The instantaneous firing rate to the 
stimulations described in D calculated from the first inter-spike interval depends on past activities. Dashed lines correspond to pre-stimulus low 
paradigm, solid lines correspond to pre-stimulus high. 
doi:10.1371/journal.pcbi.1003526.g002 



previous paragraph we calculated the m.stantaneous firing rate of 
the model when delivering a supra-threshold alternating square 
current. As expected, when a =1.0 the spike frequency shows 
adaptation only due to the membrane time constant (not shown). 
As shown experimentally and for values of a < 1 , when the input 
current makes a transition from a low input state to a high input 
state the spiking activity increases and then relaxes over time to a 
lower firing rate. We caU this downward spike adaptation. 
Similarly, when the input goes from a high to low state the firing 
rate decreases and then adapts back to a higher firing rate. We call 
this upward spike adaptation (Fig. 3H). As done experimentally, 



we fitted a single exponential to the downward or upward 
adaptation responses. The results show that the time constants of 
the upward and downward adaptation are constant across the 
multiple periods of the input signal when the signal has a fixed 
period (not shown). However, again as shown experimentally, the 
time constant of adaptation is a function of the period of the input 
(Fig. 31 and J) [37]. Our model further replicates the firing rate 
adaptation reported by the same group for neurons receiving an 
input of period 16 s (Fig. IC in [37]). In their report the authors 
calculated an average value of a ~ 0.1 5 which is well fitted by our 
model using the same value obtained by a least-squares fit 
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Figure 3. Sub-threshold and spiking fractional dynamics to oscillatory inputs. (A-B) Voltage responses, impedances, and phase angles of 
the fractional order model in response to a ZAP current. (A) A time varying sub-threshold current input (Top) and the voltage responses for three 
different values of a (Bottom). (B) Impedance and phase analysis for the simulations in A as a function of a. (C-D) Spiking response to just above 
threshold sinusoidal input (C, Top). The neuron generates an increasing number of spikes per cycle (C, Bottom). (D) The number of spikes per cycle for 
identical input as in C and varying a. (E-G) Response to supra-threshold sinusoidal input. (E) The fractional model with s( = 0.15 instantaneous firing 
rate (black) in response to a sine wave input (blue). (F) The gain of the firing rate with respect to the period length of the input (Top) shows power- 
law dynamics when plotted in log-log (Bottom). The slope of the best-fit line (red) for the log-log gain curve is -s(g,„„. (G) The phase lead of the firing 
rate in response to the same sine wave current. (H-J) Response to square periodic input. (H) In response to square wave current (Top), the fractional 
model displays upward and downward spike rate adaptation for a< 1. (I) The instantaneous firing rate shows upward and downward adaptations in 
response to changes in the period of the square input (4, 8 and 16 s). (J) The time constants of both upward and downward adaptations in (I) increase 
when the period of the alternating input current increases. (K) The spike rate adaptation of L2/3 neocortical pyramidal neurons with period 16 s 
(Fig. 1C in [37]) is fitted with the spike rate adaptation of the fractional model with a = 0.15 + 0.01 (95% confidence interval) using least-squares 
fitting. The alternating input current is switched between 3.4 and 4 nA. For E-K we used Trcr = 8 ms and t„, = 30 ms to better replicate the 
experimental data. 

doi:10.1371/journal.pcbi.1003526.g003 



(Fig. 3K). Our modeling results indicate that the fractional order of 
the derivative is less than 0.2 in order to replicate experimental 
results, which is in good agreement with the fractional exponent 
determined experimentally [37]. 

History dependent spike adaptation. In this section we 
analyze the adaptation and history dependence of the ISI. The 
spiking dynamics are affected by the long-range dependence of the 
voltage dynamics and the intrinsic membrane conductances. Both 
long-term and short-term intrinsic memory traces are observed in 
pyramidal neurons [37,38,59]. It has been recently found that the 
spiking activity of some Layer 5 pyramidal neurons in primary 



motor cortex has a long first spike latency and its ISIs decrease 
over time [38]. Using standard parameters we varied the value of a 
when stimulating the model with a step current (Fig. 4A left). We 
then fitted the plot of the ISI versus the ISI number to the data 
reported by Miller et al. (Fig. 2B in [38]). We found that when 
a<0.2 the dynamics of the ISI can be reproduced (Fig. 4A right). 
The same authors studied adaptation by applying a supra- 
threshold current for five cycles separated by an inter-stimulus 
interval (gap) [38]. We implemented this protocol in our model 
and compared the dynamics of the inter-spike intervals during 
Cycle 1 and Cycle 5. When the gap is as short as 100 ms, our 
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results show that there is spike adaptation during Cycle 1 but very 
little in Cycle 5 (Fig. 4B). In contrast, when the gap is increased to 
longer durations the spiking dynamics of Cycle 5 is very similar to 
that of Cycle 1 (Fig. 4C). We calculated the ratio of the second-to- 
last ISI in Cycles 1 and 5 to quantify spike adaptation as a function 
of the gap. Using this metric Miller et al. found that total recovery 
required a gap of 1 s (their Fig. 7C [38]); in our case the recovery 
takes about 1.5 s. However, absolute recovery of the first ISI 
between Cycles 1 and 5 took longer than 25 s (not shown). Similar 
to spike adaptation to constant input, neurons respond with a 
silent period after strong stimulation [18]. We rephcated this 
phenomenon by decreasing a to 0. 1 and injecting a square current 
with different amplitudes on top of a constant supra-threshold 
input for a fixed period of time (Fig. 4D left). After the square input 
was turned off there was a pause in the spiking activity that 
depended on the input strength (Fig. 4D right). Of course, since 
the maximum firing rate is bounded by the refractory period then 
the system reaches a limit (around 950 ms and 9 nA). This 
qualitative behavior of the fractional model is in good agreement 
with the behavior of the integrate-and-fire model with adaptation 
current (Fig. 5 in [18]), and with the firing behavior observed 
experimentally in rat subthalamic neurons (Fig. 7 in [60]). Thus, 
there is a slowly decaying memory trace that can affect the ISI 
adaptation for very long periods of time. Our analyses suggest that 
all these different neurons might be fractional differentiators. 

Spike rate adaptation to changes in the variance of noisy 
input. The adapting behavior of neurons also depends on the 
variance of the noisy input current [37,61]. We injected a noisy 
current with zero mean and varying standard deviation to the 
fractional model with a = 0.1 5 (Fig. 5A top). Our results show that 
the spiking activity of the neuron has upward and downward firing 
rate adaptation (Fig. 5 A bottom) (compare with Fig. ID in [37] 
and Fig. 3A in [61]). As in the case of the oscillatory noiseless 
input, the adaptation time constant increases directiy proportional 
to the period of the input (not shown). This property allows the 
model to follow changes in input variance as opposed to the 
classical leaky integrate-and-fire model (Fig. 5B). Changing the 
variance of the input noise can affect the spike rate of the classical 
integrate and fire model only when the mean current is very close 
to the threshold current (not shown). Thus, the fractional model 
can reproduce variance-dependent adapting behaviors such as the 
ones observed in L2/3 neocortical pyramidal neurons (Fig. 5B in 
[37]) and barrel cortex neurons (Fig. 2 in [62]). 

Spike-time reliability of the fractional integrate-and-fire 
model. Unlike constant inputs, stochastically fluctuating inputs 
can generate spike trains with highly reliable spike timing across 
repetitions [63-66] . The reliability of spike patterns is influenced 
by the mean, variance, and frequency of the stochastic input [67- 
69]. To test the reliability of the fractional model we again 
injected a constant current with added Gaussian noise. The value 
of the current was chosen for each a. so that the model produced 
low firing rate for all \'alues of a (average 14 spikes/s). A raster 
plot analysis of our results shows that the inter-trial variability 
decreases as a decreases (Fig. 6A). We quantified this reliability 
using a widely used correlation-based measure [68-71] (see 
Methods). This analysis shows that the spike reliability increases 
as the fractional exponent a decreases (Fig. 6B). The reliability of 
a neuron also increases when the variance of a time-varying 
signal embedded in a noisy signal increases [63,64,69,72]. We 
added a fixed noisy signal v](t) to the injected current and variable 
noise to analyze reliability for a = 0.1. The results show that the 
spike-time reliability depends on the variance of the embedded 
signal (Fig. 6C) [63,64,69,72]. Taken with our previous results, 
our analysis shows that when a is decreased the fractional model 



produces spike trains with strong spike time adaptation and high 
spike reliability. 

The properties of the fractional derivative that provide a 
memory trace to spiking dynamics 

The voltage memory provides non-local dynamics that affects 
the spiking activity of the cell. The voltage-memory trace decays 
over time and is dependent on the value of a. For a = 1 the process 
is identical to the classical leaky integrator, while for values of a < 1 
the past trajectory activity increasingly contributes to the present 
value of the voltage. Since the weight ^^Jv(o^^fc) is always positive, 
the sign of the voltage-memory trace depends on SVitk). With 
positive applied current the fractional model generates action 
potentials (Fig. 7A), with A V(tj^ positive until the cell fires a spike 
and is reset (Fig. 7B). However, when there is a spike and voltage is 
reset, AF(ffc) becomes negative. After the voltage escapes from the 
refractory period, the voltage-memory trace is positive until the 
next spike (Fig. 7C). As opposed to the classical leaky integrate- 
and-fire model, the memory trace accumulates over multiple 
spiking events, changes its dynamics, and thus it affects the ISI 
(Fig. 7D). 

The weight of the voltage-memory trace Wff{a,k) is determined 
by the fractional order a.. The weight Wff(ix,k) is 0 for a = 1 and it 
is always positive for o: < 1 . Fig. 7E shows the results of a 
simulation with standard parameters for 100 time steps (A^= 100). 
The X-axis corresponds to the 1-100 temporal weights at t = N. 
The y-axis corresponds to the value of each weight fFjv(a,A:). The 
increase in weight can be interpreted as the influence of the past 
state on the future state of the voltage. In a classical leaky integrator 
any past value is forgotten as a function of the time constant. In the 
fractional leaky integrator all the past values could continue to 
influence the future behavior of the system, particularly, for low 
values of a. Fig. 7F illustrates this point by showing the dynamics of 
the weight of the initial condition (Wf^iafi)) as a function of a. The 
other important term that comes from the fractional derivative is the 
fractional coefficient K(a,dt) = (dtfr(2 — a) which is a weight 
factor for the Markov process (Eq. 6). When a becomes smaller, this 
function grows rapidly and it makes the effect of the input current 
on the voltage dynamics stronger (Fig. 7G). It is the combination of 
the weighted Markov process and the opposite effect of the memory 
trace that contribute to the long term spiking dynamics of the 
fractional model. 

Most of our results suggest that the value of a has to be low in 
order to reproduce the spike timing adaptation observed experi- 
mentally. Although the fractional model has a continuous 
dependence on a the power-law dynamics cause the effects to be 
nonlinear. For oc close to 1 the effects of the Markov term weighted 
by the gamma function dominate the dynamics. It is only when a 
decreases that the voltage-memory trace can slow down the time 
evolution of the voltage. This is illustrated by plotting the value of 
the Markov term versus the memory trace for simulations in which 
we apply a step current (Fig. 8). When 0!= 1 the memory trace is 
zero and tlu' \ oltage only moves along the Markov term axis. As the 
value of a decreases the voltage trajectory is deflected, taking longer 
to depolarize. The power-law dependency results that when the 
value of a < 0.2 in that the memory trace dominates in the initial 
moments of the depolarization and slows down the dynamics (Fig. 8). 

Comparison to other models of fractional integrate-and- 
fire 

When the input current is constant and o: = 1 , the sub-threshold 
voltage dynamics (Eq. 1) have an analytic solution which is the 
same as the solution of the classical integrate-and-fire model [73]. 
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Figure 4. Inter-spike interval adaptation and history dependence. (A) Left: The spiking activity of the fractional model with a step current of 
4 nA. Right: The inter-spike interval (ISI) curve of Layer 5 pyramidal neurons in primary motor cortex (Fig. 2B in [38]) is fitted with the ISI curve of the 
fractional model with a = 0.19 + 0.02 (95% confidence interval with a least-squares fit). The first 7 ISIs of the model are removed for best fit. (B-C) 
IVIodeling the intrinsic memory of adapting pyramidal cells. (B) Left: The voltage trace of the model in response to a step current separated by 0.1 s 
inter-stimulus interval. Right: The ISIs of Cycle 1 and 5 as a function of ISI number. (C) The same as panel B, but with longer inter-stimulus interval, 2 s. 
For A-C a = 0.19. (D) Memory induced pauses of the model with a = 0.1 depend on the magnitude of the current pulse. Left: Voltage traces with 
shorter and longer pauses in response to 1 nA and 4 nA current pulses, respectively. Right: The pause of the spiking activity increases as a function of 
the magnitude of the current pulse. For all we used trcf = 8 ms and t,„ = 30 ms to better replicate the experimental data. 
doi:10.1371/journal.pcbi.1003526.g004 



Similarly, for constant input current Langlands et al. [74] derived 
the analytic solution of the sub- threshold voltage for a < 1 from the 
fractional cable equation. In that model the integration of the 
memory trace is restarted after every spike, thus wiping out the 
memory trace (see Methods). Such a system is capable of 
reproducing the delay to first-spike to constant input but then 
produces regular spiking (Fig. 9A). Our fractional model replicates 
the same result when we reset the memory trace after every spike 
(Fig. 9B). However, our model greatly differs from the analytical 
solution when taking into account the cumulative effect of the 
memory trace across multiple spiking cycles (Fig. 9C-D). Although 
both the analytic and full model with memory reset can capture 
short term memory, they do not produce spike adaptation. Hence, 
the full fractional model without any memory reset captures the 
multi-scale processes that spans the spiking activity of neurons. 



Discussion 

Spike timing adaptation is a widespread phenomenon through- 
out the nervous system [37,39,75]. In particular, neocortical 
pyramidal cells produce spike adaptation with multiple timescale 
dynamics [5,37,59]. Our model is capable of reproducing multiple 
sub-threshold and spike timing adaptation properties reported by 
different groups and with different experimental conditions. The 
conclusion from fitting our model to experimental results is that 
a < 0.2. This indicates that the order of the fractional derivative has 
to be very low for the memory trace to overcome the classical 
contribution of the leaky integrator. Furthermore, the fractional 
model is capable of producing spike trains with high adaption and 
reliability. Our work provides a framework to study spike adaptation 
as part of power-law dynamics and the techniques used here can be 



PLOS Computational Biology | www.ploscompbiol.org 



8 



March 2014 | Volume 10 | Issue 3 | el 003526 



Fractional Integrate-and-Fire Model 




Time (s) Time (s) 

Figure 5. Firing rate adaptation to cKianges in input variance. (A) Top: Noisy input current with two standard deviations: /i„j =6+^(^(0 nA, 

A = ^ or 2 nA. Bottom: The firing rate to noisy input current calculated from 100 trials, 2 = 0.15. (B) Top: Time varying noisy input current. /i„j = 5+ 
Ai;{t) nA, ^ = 1, 4, 2, 1, 2, 1, 4 and 2 nA, consecutively, and the noise c(t) is filtered with an alpha function /(() = (e'^'/^' with t = 0.2 ms. Bottom: 
Instantaneous firing rate in response to the input current for o!=1.00 (blue) and s( = 0.15 (black). For both Tref = 8 ms and t„, = 30 ms. 
doi:1 0.1 371 /journal.pcbi.1 003526.g005 



applied experimentally to determine if a neuron is following power- 
law adaptation from the sub-threshold to firing rate regimes. 

Using fractional dynamics to study non-local interactions 

The fractional model can produce different degrees of adapting 
electrical activities by modifying the fractional exponent a. A 
fractional order derivative captures the long-range correlations of 
the a system models that results in non-local dynamics. Fractional 
differential equations have been used in biological systems 
to capture the long-term memory effects of the dynamics 
[33,74,76-81]. For example, fractional order derivatives have 
been observed in the vestibular-ocular system [27,82] and in the 
gating dynamics of ion channels [44,83]. Power law statistical 
distributions of a neuronal response also exhibit fractional order 



dynamics [37]. The voltage dynamics in the fractional order 
model depend on both the Markov term (immediate past) and the 
voltage-memory trace that integrates all past voltage values. The 
behavior of the voltage-memory trace is similar to the behavior of 
the adaptive filter in the work of Pozzorini et al. [59], although in 
their work this filter is described as the sum of the spike-triggered 
current and a moving threshold. The voltage-memory trace also 
corresponds to the adaptation integral used in other works [18]. In 
this context our fractional order leaky integrate-and-frre model is a 
unified mathematical and computational framework that can be 
used to describe power-law dynamics and long-range correlations 
in neuronal activity. The fractional derivative can capture 
relationships between the distribution of conductances that can 
be complicated to model using explicit techniques. 
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Figure 6. Spil<e-time reliability increases as the fractional exponent a decreases. (A) Raster plots of the response of the fractional model to 
a noisy input under three different values of a. (B) Spike-time reliability of the fractional model increases as the fractional exponent a. decreases. (C) 
Reliability increases when the standard deviation of an embedded fixed signal increases. See Text and Methods for details. 
doi:1 0.1 371/journal.pcbi.1 003526.g006 
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Figure 7. The properties of tlie voltage-memory trace. (A-D)The changing response of the memory trace across multiple spikes, a = 0.2. (A) 
Voltage trace of the fractional model stimulated with /j„j = 3 nA. Spikes have been clipped to emphasize the sub-threshold dynamics. (B) AK for the 
data in A. (C) The memory trace for the data in A. (D) Overlapped memory traces for different inter-spike intervals during the same simulation. (E-G) 
The dynamics of the weight of the voltage-memory trace fV^i'^,!^) and the fractional coefficient depend on the fractional exponent a. (E) When a 
decreases the weights increase. (F) The value of the weight IVn{3.,0} as a function of a. and time. (G) The fractional coefficient of the Markov process 
K{a.,dl) increases when i is decreased. 
doi:10.1371/journal.pcbi.1003526.g007 



The biophysical interpretation of the fractional order of a 
neuron 

In biophysical systems, the voltage-memory trace might 
represent spike-triggered mechanisms that cause adaptation. For 
example, the voltage-memory trace might represent the slowly 
inactivating potassium-like current that induces the upward spike 
adaptation shown in Layer V pyramidal neurons of primary motor 
cortex [38], suggesting that these neurons are fractional differ- 
entiators. Alternatively, the voltage-memory trace can correspond 
to other adaptation currents such as calcium-activated after- 
hyperpolarization currents [5], slow sodium-channel inactivation 
currents [84,85], or a combination of several adaptation currents. 
Many studies have used models with slow adaptation currents and 
exponential functions to analyze spike time adaptation [5,54,84]. 



a =1.0 




-70 L. f I I I I I I — I 
0 0.1 0.2 0.3 
Memory trace (mV) 

Figure 8. The memory trace dominates tKie fractional dynamics 
for low values of a. Markov term versus memory trace as a function of 
a. The fractional model was stimulated with constant current (0.3 nA) 
for 5 seconds. For d;= 1 the memory trace is zero and the voltage only 
moves along the y-axis. 
doi:1 0.1 371/journal.pcbi.1 003526.g008 



Some of these models can produce similar properties found in 
fractional dynamics. For example, the Hodgkin-Huxley model 
with slow after hyper-polarization currents can produce multiple 
timescale adaptation processes [37]. The Generalized leaky 
integrate-and-fire model with an adaptive filter (GLIF-^) produces 
spike adaptation with power-law dynamics [59]. Both the power- 
law dynamics and history-dependent properties of the GLIF-(^ 
model correspond to that of the fractional leaky integrate-and-fire 
model. However, the fractional model provides a general way by 
simplifying the complicated details shown in other models. The 
fractional model exhibits spike adaptation with power-law 
dynamics by integrating all the past voltage values without any 
additional adaptation currents. Power law functions generalize the 
mechanism underlying exponential processes and are better 
alternatives to describe scale invariant spike adaptation 
[17,18,42,59,86,87]. The fractional model shows new directions 
for studying spike adaptation using fractional derivatives and 
power-law dynamics instead of classical derivatives and exponen- 
tial functions. 

The key parameter in our fractional model is a. Experimental 
results have suggested that a. can be as small as 0.15 for 
neocortical pyramidal neurons [37]. Our fitting of these data also 
resulted in a value of a « 0.1 5. Fits to the response of Layer 5 
pyramidal motor cortex neurons resulted in a value of a « 0.1 9 
[38]. Thus, the experimental results and our modeling analysis 
suggest that the biophysically important values of a are when 
a<0.2. However, it is clear that if a is much closer to 0, the 
system takes a longer time to generate spikes or never fires spikes 
at all, depending on the magnitude of the stimulus, so the feasible 
range of the fractional exponent might be between 0.05 and 0.2. 
The value of a might correspond to the type, function, or location 
of specific neurons [38]. These regional differences are not 
exclusive to the cortex, for example, Purkinje cells in Lobule X 
and Lobules III - V show different degrees of spike adaptation 
[75]. Thus, different values of a can be used to map the general 
voltage and spike time adaptation properties of neurons 
throughout the brain. 
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Figure 9. The fractional model and its analytic solution with memory reset show no spike adaptation. (A) The spike train produced by 
the analytic solution with memory reset displays regular spiking. (B) The spike train produced with the fractional model with memory reset also 
displays regular spiking. (C) The spike train produced by the full fractional model without any memory reset display spike adaptation. (D). The firing 



rate curves of the analytic solution, fractional model with memory reset and full fractional model. For all panels !)! = 0.1, /jn 

Tref = 8 ms. 

doi:1 0.1 371/journal.pcbi.1 003526.g009 



nA, rm = 30 ms and 



Determining if a neuron is a fractional differentiator 

Previous work ha.s provided an experimental foundation to 
determine if a spiking neuron is a fractional differentiator [37]. 
Our framework provides a more general methodology to 
determine power-law neuronal dynamics from the sub-threshold 
to the spiking regime. At the sub-threshold level there are several 
measurements that can indicate that the membrane of the neuron 
follows a power-law. For example: 

• The flattening of the impedance and phase response to the 
ZAP current. 

• A power-law behavior can be measured by applying a low 
depolarizing current and plotting Log(voltage/time) versus 
Log(time). A straight oblique line indicates power-law 
behavior. 

In the same experiment a series of protocols can also be applied 
to determine if the spiking activity follows power-law dynamics. 
Some of these measurements are straightforward, others require 
longer recordings. For instance: 

• The ISI histogram plotted as Log(Counts) versus Log(ISI) 
follows a power-law. 

• The instantaneous firing rate response to step currents depends 
on the value of the current before the stimulation. 

• Using a sinusoidal or square oscillatory input with varying 
period, fit a time constant to the adapting firing rate. If the 
time constant depends on the stimulus this suggests a memory 
trace. 

• The gain of the spiking neuron follows a power-law. 

• The ISI adaptation depends on the inter-cycle (gap) time and 
the pause length of the neuron depends on the strength of the 
previous stimulus. 



• The neuron shows firing rate adaptation to changes in 
variance with frxed mean. 

The computational importance of power-law spike time 
adaptation 

Although power-laws are found at multiple scales of biological 
organization, their function and importance are still debated 
[88,89]. In our work we propose that the membrane voltage of 
neurons can follow a power-law due to the emergent property of 
the combination of multiple active conductances. The value of the 
fractional derivative can be mapped to spike time adaptation 
dynamics taking place in multiple cell types across the brain. 
Computationally, a low value of a. results in spiking dynamics that 
are at the same time highly adaptable and reliable. Thus, neurons 
following power-law adaptation could have a large operational 
range while providing the reliability to filter out noisy signals and 
increase information capacity. The lack of a sub-threshold 
resonance frequency allows the neuron to fdter signals homoge- 
neously over a wide range of frequencies. In such a case, the 
fractional leaky integrate-and-fire model provides the basis to 
study the computational capacities and information processing 
properties of neurons showing high degree of spike time 
adaptation. 

Methods 

Implementing the fractional leaky integrate-and-fire 
model 

The equations were coded and implemented using our recendy 
developed fractional integration toolbox [47] and the simulation 
software package MATLAB [90] . The toolbox can be downloaded 
at www.utsa.edu/SantamariaLab. The parameters for all simula- 
tions were frxed and are described in Table 1. 
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Comparing models to experiments 

In order to compare to experiments we extracted the data 
from the referenced material using WebPlotDigitizer (http:// 
arohatgi.info/WebPlotDigitizer/). Then we imported the data 
points into Matiab. We then ran simulations varying the value of 
a, usually between [0.5, 1.0] at 0.1 steps. We minimized the 
mean squared error between the data and the simulations, 

E = (l /N)^Jj2 (data — simulation)^, where N is the number 
of points. We determined the 95% confidence intervals by 
then varying a around this minimum value of E and 
calculated when it changed for an amount larger than 5% in 
either direction. 

Modeling oscillatory input 

In order to get only sub-threshold oscillations we used 

7i„j=0.3sin(27c/0 where /= 10o(^ ^ _^_^(^,/ i5oo) 

sigmoidal frequency function of time that \-aries from 0 Hz to 
100 Hz in 10 seconds. The impedance as a function of frequency 
is defined as Z(f)= V(F)/ I(f) = R(f)+jX(f), where V is the 
membrane voltage, I is the input current, R is the resistance and X 
is the reactance. The absolute value of Z can be calculated using a 
fast Fourier transform \Z\ = \FFT( V)/FFT{r)\ and the phase as 
phase =a.Tctan(X{f)/R(f)) [57]. 

Reliability quantification 

The simulations to study spike time reliability were generated by 
injecting a current Iiaj=I + A^(t) nA, where ^(t) is a Gaussian 
white noise with zero mean and standard deviation A = 0.03 nA. 
The stochastic input t^(/) is filtered with an alpha function 
f{t) = te^'/^ with time constant t = 3 ms. The spike trains 
were obtained from A'^ = 50 trials, and the trail-to-trail variability 
of those N different responses were caused by the noise i,{t) while 7 
and A were fixed [68,71]. In order to avoid initial condition 
effects we analyzed the spike trains of the last 5 s from 10 s 
simulations. 

The reliability measurements were computing using a correla- 
tion-based measure [68-71]. In brief, the spike trains obtained 
from N trails were smoothed with a Gaussian filter of width 3(7^) 
and then pairwise correlated. The correlation-based measure 
reliability R is defined as 



R-- 



N{N-\) 



EE 

i=l y=,+l 



\t\\Sj\ 



(10) 



where N is the number of trials and the vectors Sj{i=\,...,N) 
are the filtered spike trains, filtered using Oc = 3 ms. The values 
of R range from 0 (lowest reliability) to 1 (highest reliability), 
and the reliability R was computed for a in the range [0.1, 
1.0]. 

For the quantification of the rehability in the coding of 

an embedded signal we injected the following current 
hn)=I ''r A^(i) + afr](t) nA, where 7 = 5 nA is the mean current 
which generates very low firing rate, A = \ nA is the standard 
deviation of the intrinsic noise, and oy is the standard devia- 
tion of the embedded signal also generated by a Gaussian 
noise and varies from 0 to 6 nA. The intrinsic and embedded 
signals are filtered with an alpha function f(t)= —te^'l'' where 
T = 3 ms. 



The comparison of the analytical and full fractional 
model 

We compared our fractional model to a previously developed 
analytical model of a fractional leaky integrate-and-fire [74]. 
Briefly, this analytical model is obtained with the following steps. 
Equation 1 can be converted to 



^ = ^((K-K.)+^), (11) 

af" -Cm \ gLj 

where = — is the membrane time constant. By applying . _ 
on both sides we obtained 



dt T„ dt^ "V Sl 



(12) 



Equation 12 is solved using the Fourier-Laplace transform (for 
details see [36,74]) and the solution is given by 



nty- 



gL 



V{to)-Vj^-^f^E,[-^^-^], (13) 



where E^ is the Mittag-Lefiler function [36], and for small times 
this function is approximated as 



eA - 



{t-t<,r 



iexp 



(14) 



For the simulation of this model presented in our work we used 
the fuU Mittag-Lefiler function instead of this approximation. In 
Eq. 13 the term with the Mittag-Lefiler function (right side and 
right term) represents the memory trace. As in the classical 
integrate and fire model when the voltage reaches Ftj, a spike is 
generated and V is reset to Freset for a refractory period Tref . The 
subthreshold voltage is integrated using Eq. 1 3 with initial voltage 
V{to) = I^reset and new initial time to. During each integration 
cycle the Mittag-Lefiler function restarts from 0 since the initial 
time ?o reset to a new value. We call this memor\' reset. Thus, this 
model wipes out the memory trace after every spike, in contrast to 
our model that integrates the entire voltage trace. As a result, the 
inter-spike intervals of the spike train of the analytic solution are 
equal. If the approximation (Eq. 14) is used to simulate the voltage, 
the firing rate can also be approximated analytically by combining 
Eq. 13 and Eq. 14 (see also [74]). Let 7/ be the time when the 
voltage takes to increase from Freset to Fth and to fire. The time Tf 
is given by [74] 



■Vl- 



gL 



V,h-Vl- 



gL 



(15) 
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Using the above the firing rate is approximated by 

firing rate « . (16) 

In the Results section we compare this analytical model with our 
model with memory reset (re-starting the memory trace after every 
spike) and with the full model (integrating the memory trace from 
the beginning of the simulation). 
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